Mapping

Basic map

Add elements to map

plot(wrld_simpl,xlim=c(115,128) ,ylim=c(19.5,27.5),col='#D2B48C',bg='lightblue') # TW map
coords <- matrix(c(121.537290,120.265541, 25.021335, 22.626524),ncol=2) # NTU and SYS univ. 
coords <- coordinates(coords) # assign values as spatial coordinates
spoints <- SpatialPoints(coords) # create SpatialPoints
df <- data.frame(location=c("NTU","SYS")) # create a dataframe
spointsdf <- SpatialPointsDataFrame(spoints,df) # create a SpatialPointsDataFrame
plot(spointsdf,add=T,col=c('white','red'),pch=19,cex=0.5) # plot it on our map
text(121,24, 'TAIWAN', cex=1)

• Coordinates of these two university are not in spatial coordinate format, which means it may run off while you change projection style of your map

Fix position with coordinates

plot(wrld_simpl,xlim=c(-130,-60),ylim=c(45,80),col='#D2B48C',bg='lightblue')
coords <- matrix(c(-110,-102,-102,-110,-110,60,60,49,49,60),ncol=2)
l <- Line(coords)
ls <- Lines(list(l),ID="1")
sls <- SpatialLines(list(ls))
df <- data.frame(province="Saskatchewan") #create an empty dataframe to put lines into it.
sldf <- SpatialLinesDataFrame(sls,df)
plot(sldf,add=T,col='#3d2402', cex=2)
text(-114, 55, 'Saskatchewan', srt=90, cex=0.5)
text(-114, 63, 'CANADA', cex=1)

generate map from GADM data

library(raster)
TWN1 <- getData('GADM', country="TWN", level=0) # data Taiwan

## Warning in getData("GADM", country = "TWN", level = 0): getData will be removed in a future version of raster
## . Please use the geodata package instead

JPN <- getData('GADM', country="JPN", level=0) # data Japan

## Warning in getData("GADM", country = "JPN", level = 0): getData will be removed in a future version of raster
## . Please use the geodata package instead

class(TWN1) # those datasets are SpatialPolygonsDataFrame

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

par(mfrow = c(1, 2))
plot(TWN1,axes=T,bg=colors()[431],col='grey')
plot(JPN,axes=T,bg=colors()[431],col='grey')

Zoom on map

plot (TWN1, axes=T, xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431],col='gray')

plot (TWN1, axes=T, xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431],col='grey') 

TWN2 <- getData('GADM', country="TWN", level=1)

## Warning in getData("GADM", country = "TWN", level = 1): getData will be removed in a future version of raster
## . Please use the geodata package instead

TWN2$NAME_1

## [1] "Fujian"     "Kaohsiung"  "New Taipei" "Taichung"   "Tainan"    
## [6] "Taipei"     "Taiwan"

plot(TWN1,col="grey",xlim=c(119,122.5), ylim=c(21.5,25.5), bg=colors()[431], axes=T)
KAO <- TWN2[TWN2$NAME_1=="Kaohsiung",]
plot(KAO,col="grey 33",add=TRUE)

# base map
plot(TWN1,col="grey",xlim=c(121,122), ylim=c(24,25.5), bg=colors()[431], axes=T)
# adding  spatial polygones 
TAI <- TWN2[TWN2$NAME_1=="Taipei" | TWN2$NAME_1=="New Taipei" ,]
plot(TAI,col="black",add=TRUE)
# adding spatial points 
coords <- matrix(cbind(lon=c(121.2,121.55,121.8),lat=c(25,25.19,24.5)),ncol=2)
coords <- coordinates(coords)
spoints <- SpatialPoints(coords)
df <- data.frame(location=c("City 1","City 2","City 3"),pop=c(138644,390095,34562))
spointsdf <- SpatialPointsDataFrame(spoints,df)
scalefactor <- sqrt(spointsdf$pop)/sqrt(max(spointsdf$pop))
plot(spointsdf,add=TRUE,col='white',pch=1,cex=scalefactor*3,lwd=2) 
# adding a location of NTU (not spatial point here)
points(121.537290,25.021335, type="p", pch=18, col='white', cex=1.5)
# adding text
text(121.53,24.921335,"NTU", col='white', font=2)
# adding scale
maps::map.scale(x=120, y=25.4)

ggplot2

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ dplyr::select()  masks raster::select()

library(rnaturalearth)
library(rnaturalearthdata)
library(sf)

## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

theme_set(theme_bw()) 
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)

## [1] "sf"         "data.frame"

ggplot(data = world) +
  geom_sf() +
    coord_sf(expand = FALSE) +
    xlab("Longitude") + ylab("Latitude") +
    ggtitle("World map", subtitle = paste0("(", length(unique(world$name)), " countries)")) +
  geom_sf(color = "black", fill = "lightgreen")

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

dev.off()

## null device 
##           1

GBIF data

library(rgbif)
library(mapr)

gbif.res <- occ_search(scientificName = "Urocissa caerulea", limit = 1200)

map_ggplot(gbif.res) +
  coord_sf(xlim = c(120, 123), ylim = c(21, 26), expand = FALSE)

Bathymetric map with marmap

library(marmap)

## Registered S3 methods overwritten by 'adehabitatMA':
##   method                       from
##   print.SpatialPixelsDataFrame sp  
##   print.SpatialPixels          sp

## 
## 載入套件：'marmap'

## 下列物件被遮斷自 'package:raster':
## 
##     as.raster

## 下列物件被遮斷自 'package:grDevices':
## 
##     as.raster

# query
TW.bathy <- getNOAA.bathy(lon1=118,lon2=124, lat1=21,lat2=26,resolution=1) # don't put too wide / resolution: 1 

## Querying NOAA database ...

## This may take seconds to minutes, depending on grid size

## Building bathy matrix ...

# define palette
blues <- colorRampPalette(c("darkblue", "cyan"))
greys <- colorRampPalette(c(grey(0.4),grey(0.99)))
# make the plot, step: difine line where to exist (e.g. draw lines per 2000 meter between -6000 to -1000 meters)
plot.bathy(TW.bathy,
     image=T,
     deepest.isobath=c(-6000,-120,-30,0),
     shallowest.isobath=c(-1000,-60,0,0),
     step=c(2000,60,30,0),
     lwd=c(0.3,1,1,2),
     lty=c(1,1,1,1),
     col=c("grey","black","black","black"), 
     drawlabels=c(T,T,T,F),
     bpal = list(c(0,max(TW.bathy),greys(100)),c(min(TW.bathy),0,blues(100))),
     land=T, xaxs="i"
     )

Leaflet

library(leaflet)
library(rgdal)

## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/ASUS/AppData/Local/R/win-library/4.2/sf/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/ASUS/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

taiwan <- readOGR('D:/R_WD/Git linked/2022RisFUN-Mo/Mapping/mapdata202209220943/COUNTY_MOI_1090820.shp', use_iconv=TRUE, encoding='UTF-8')

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum Taiwan_Datum_1997 in Proj4 definition: +proj=longlat
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\R_WD\Git linked\2022RisFUN-Mo\Mapping\mapdata202209220943\COUNTY_MOI_1090820.shp", layer: "COUNTY_MOI_1090820"
## with 22 features
## It has 4 fields

FRE <- paste(sep = "<br/>",
  "<b><a href='https://www.dipintothereef.com/'>FRELAb TAIWAN</a></b>",
  "Functional Reef Ecology Lab",
  "Institute of Oceanography, NTU")

# Change style by addtiles
leaflet(taiwan) %>%
  addPolygons(weight=0.5) %>%
  addTiles(group="Kort") %>%
  addPopups(121.53725, 25.021252, FRE, options = popupOptions(closeButton = FALSE))